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Modelling the Bioclimatic Niche and Distribution of the Steppe Mouse, Mus spicilegus (Rodentia, 
Muridae), in Ukraine. Tytar, V. M., Kozinenko, I. I., Mezhzherin, S. V. — The Steppe mouse, Mus 
spicilegus, is endemic to Europe and found to be expanding its home range in recent years. In Ukraine there 
are indications a north- and eastwards expansion and/or reestablishment of M. spicilegus. We suggest that 
climatic conditions may be the primary factors that foster or limit the range expansion of M. spicilegus 
in Eastern Europe. Our objective was to complement the knowledge about the distribution of the species 
with an estimation of the potential distribution of the species in Ukraine using known occurrence sites 
(in Ukraine and neighbouring areas) and environmental variables in an ecological niche modelling 
algorithm. After accounting for sampling bias and spatial autocorrelation, we retained 73 occurrence 
records. The algorithm used in this paper, Maxent (Phillips et al., 2006), is a machine learning algorithm 
and only needs presence data, besides the environmental layers. Using this approach, we have highlighted 
the importance and significance of a number of bioclimatic variables, particularly those characterizing 
wintering conditions, under which higher mean temperatures enhance habitat suitability, whereas 
increased precipitation leads to an opposite effect. The broadly northwards shift of the home range of the 
species in Ukraine could generally be due to the increasing (since the 1980s) mean temperature of the 
winter season. We expect this expansion process will continue together with the changing climate and 
new records of locations of the species may be used for monitoring such change. 
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Introduction 


The Hillock, or Mound-building, or Steppe mouse, Mus spicilegus, is endemic to Europe, found from 
Lake Neusiedl on the Austro-Hungarian border expanding through Slovakia, Hungary, Serbia, Montenegro, 
Bulgaria, Moldova, and Ukraine, reaching as far as the Rostov-on-Don region in the utmost south-west of 
Russia (Coroiu et al., 2016). It is unique among mice in its habit of building earthern mounds constructed 
between about mid-August and mid-November. The construction is done by juveniles, three to four weeks 
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old, and the mound is their winter home, with adults occasionally cohabiting (Sokolov et al., 2008), however 
by spring adults have died, leaving only the off-spring to found the following generation (Sokolov et al., 1990). 
The mouse occurs in a variety of habitats including grassland, steppe, cultivated land, orchards, clearings and 
woodland borders (Coroiu et al., 2016). 

The Steppe mouse is a common species with a wide range, and the International Union for Conservation 
of Nature (IUCN) has rated its conservation status as being of “least concern” (LC) while noting that the 
intensification of agriculture and destruction of grassy steppe may be a future threat. However, in Slovakia the 
range appears to be expanding (Coroiu et al., 2016). The same seems to be occurring in Ukraine and Russia. 
Recent records are indicating a north- and eastwards expansion of the species’ range and/or reestablishment 
of M. spicilegus in places from which the species had disappeared since the first half of the XXth century. For 
instance, a thorough search of the species in the Ukrainian steppe reserve of Askania-Nova found no trace of 
M. spicilegus until 1999 (Polishchuk, 2012), despite having records published early as 1928 (Brauner, 1928), 
and the same relates to the Rostov-on-Don region of Russia, where the species reappeared in 1999-2004 in 
seven locations (Lipkovich, 2005). On the other hand, M. spicilegus had never previously been recorded from 
the Kaniv Nature Reserve (49.74444° N, 31.45583° E), a location nearby the northern boundary of the species’ 
home range, so the findings made here in 2000-2004 (Ruzhilenko, 2005) could be evidence of range expansion, 
which, in the opinion of the author, are due to climate change. 

According to I. V. Zagorodnuik, the range of the species in Ukraine is split into parts by at least four 
large rivers, the Prut, Dnister, Southern Boug and Dnipro, as well as by an extensive number of smaller ones, 
whose sources are located considerably further to the north beyond the contemporary range of M. spicilegus. 
Despite this, in all cases the northernmost findings of the species are located on the same line; see figure 2 in 
Zagorodnuik, Berezovsky (1994). Because winter migrations of the animals are hardly possible, obviously there 
are specific external factors that determine the character of the species’ geographical distribution. Obvious 
biogeographic barriers, as noted above, are absent. There are also no obvious differences in the habitat (moisture, 
biotope spectrum and their feeding capacity) within the distribution area of the species in more northern areas. 
Therefore, as I. V. Zagorodnuik concludes, there is a need to search for other reasons, among which there may 
be wintering conditions for the species, primarily the temperature and depth of soil freezing (Zagorodnuik, 
Berezovsky, 1994). Earlier an assumption had been made that Steppe mouse population reductions are largely 
due to adult deaths occurring during the winter (Berry, 1981). Later on the main factors responsible for mouse 
deaths were concretized as low temperature and frequent thaws, leading to drenched food and nest chambers 
(Muntyanu, 1990). 

Under the current knowledge, we in our turn suggest that climatic conditions may be the primary 
factors that foster or limit the range expansion of M. spicilegus in Eastern Europe. Our objective was to 
complement the knowledge about the distribution of the species with an estimation of the potential 
distribution of the species in Ukraine using known occurrence sites (in Ukraine and neighbouring areas) 
and environmental variables in an ecological niche modeling algorithm. Ecological niche models (ENMs) 
have been widely used for predicting potential distributions of species and exploring their habitat suitability 
requirements (Franklin, 2009; Miller, 2010). This methodology can be particularly useful to predict the 
potential distribution of M. spicilegus because climatic conditions are likely to play an important role in 
limiting the distribution of this species. 

Traditionally, determining such driving factors would require laborious field measurements of the key 
environmental variables in natural populations. However, the advent of GIS and the increased availability of 
global environmental data in recent years have favoured the proliferation of diverse kinds of ENMs intended 
to answer a wide range of applied ecological questions (Peterson et al., 2011). Because these models seek to 
identify the features that characterize a species’ known distribution (their “bioclimatic envelope”), ENM can 
provide basic quantitative information about species’ apparent habitat preferences (Nakazato et al., 2010). 


Material and methods 


Occurrences from the known range of the species in Ukraine, Moldova, the Rostov-on-Don and Belgorod 
regions of Russia were obtained from the literature (Lyalyukhina et al., 1989; Zagorodnuik, Berezovsky, 1994; 
Kondratenko, 1998; Ruzhilenko, 2005; Lipkovich, 2005; Tsvelykh, 2009; Smirnov, 2009; Tokarsky et al., 2011; 
Polishchuk, 2012; Evstafiev, 2015; Partolin, 2016; etc.). Field trips undertaken in 2017 to Cherkasy Region 
yielded additional, previously unknown records of the species. In total our dataset contained 269 georeferenced 
occurrences (fig. 1); occurrences that lacked latitude and longitude coordinates were georeferenced using 
Google Maps (http://maps.google.com.ua). 

These occurrence points varied in spatial density due to variable sampling intensity over geography. 
As a result, and to avoid overemphasizing heavily on sampled area, we selected points for model calibration 
using a subsampling regime to reduce sampling bias and spatial autocorrelation, which would produce models 
of lower rather than higher quality (Beck et al., 2013). Following M. A. Nuñez and K. A. Medley (Nuñez, 
Medley, 2011), we generated models using all available occurrence points and measured spatial autocorrelation 
among model pseudo-residuals (1 — probability of occurrence generated by model) by calculating Moran’s I at 
multiple distance classes using the SAM v4.0 software (Rangel et al., 2006). Significance was determined using 
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Fig. 1. Occurrences of Mus spicilegus in Ukraine and neighbouring areas used for creating the ENM. [Data 
collected before (triangles) and after (circles) 1990.] 


permutation tests (n = 999). Moran’s I is a widely used measure of spatial autocorrelation, ranging from 0 to 1, 
with values > 0.3 considered relatively large (Lichstein et al., 2002). A minimum distance of 53.6 km, at which 
Moran’s I < 0.3, p = 0.001, was detected. 

Next we used the spThin package in R (Aiello-Lammens et al., 2015) to subsample our dataset such that all 
occurrence records were separated by this minimum distance. The procedure greatly reduced sampling bias and 
spatial autocorrelation, resulting in evenly distributed occurrence points across space. In the end, we retained 
73 occurrence records after thinning. This number of occurrence records is considered more than sufficient to 
generate robust species distribution models (Hernandez et al., 2006). 

The algorithm used in this paper, Maxent (Phillips et al., 2006), has proven good performance and 
accuracy for such studies (Elith et al., 2011). Maxent (Version 3.3.3k) is a machine learning algorithm. The 
main advantage of applying Maxent to the modeling of geographical species distributions in comparison with 
other methods is that it only needs presence data, besides the environmental layers. Pseudo-absence points 
(used instead of true absences) were randomly generated (using the default option) within a bounding box 
encompassing M. spicilegus presence points. 

As part of its output, Maxent ranks the environmental layers used to train the ENM based on their 
relative importance in model formulation. To do so, it employs two metrics: (1) percent contribution, and 
(2) permutation importance. The latter is considered to provide a much more accurate ranking than percent 
contribution (Searcy, Shaffer, 2016). In the permutation importance option the contribution for each 
variable can be determined by randomly permuting the values of that variable and measuring the resulting 
model performance. Values are normalized to give percentages; we considered percentages exceeding a 10 % 
threshold. The second option to determine the importance of environmental variables uses a jackknife test 
and the regularized gain change during each iteration of the training algorithm. The environmental variable 
with the highest gain is considered to have the most useful information by itself, whereas the variable causing 
the largest decrease in the model’s gain contains the most information not found in the other environmental 
variables. 

Importantly, Maxent also allows the construction of response curves to illustrate the effect of selected 
variables on habitat suitability. These response curves consist of the specific environmental variable as the x-axis 
and, on the y-axis, the predicted probability of suitable conditions as defined by the logistic output. Upward 
trends for variables indicate a positive relationship; downward movements represent a negative relationship 
(Baldwin, 2009). 

The Maxent general-purpose machine learning technique is prone to over-fitting, therefore uses 
“regularization” (a beta parameter) to avoid over-fitting data. Higher values of beta increase the “smoothness” 


474 V. M. Tytar, I. I. Kozinenko, S. V. Mezhzherin 


of species’ responses to the environment. There is also a choice of an expanded set of transformations of the 
original covariates (termed features), which are helpful in making predictions about the outcome. MaxEnt fits 
models using linear, product, quadratic, hinge and threshold functions, therefore enhancing the flexibility for 
modelling non-linear species-habitat responses. Because MaxEnt’s default regularization parameter and feature 
classes have profound impacts on model performance (Merow et al., 2013), we used ENMeval (Muscarella et al., 
2014) to build a series of models with all possible combinations of these parameters. We selected a model with 
a combination of feature class and regularization multiplier that provided the best trade-off between model 
goodness of fit and complexity using the corrected Akaike information criterion, AICc (Akaike, 1974; Warren, 
Seifert, 2011). Replicated runs (n = 25) of bootstrap type were completed to build the model. 

We used the 10th percentile training presence logistic threshold value to generate binary maps (Liu et al., 
2005). This threshold value provides a better ecologically significant result when compared with more restricted 
threshold values (Phillips, Dudik, 2008). Based on the probability value, we divided the habitat areas into three 
classes: unsuitable and/or marginal area (below the threshold value), median area (threshold value — 0.6), and 
core area (0.6-1.0). 

Maps of habitat suitability in the ASCII format were processed and visualized in SAGA GIS (Conrad et al., 
2015). Statistical data was analyzed using the PAST software package (Hammer et al., 2001). 


Environmental variables 

For bioclimatic modeling we used the CliMond archive, a set of free climate data products (Kriticos et al., 
2012). The climate data includes interpolated surfaces at 10’ and 30’ for recent historical climate and relevant 
future climate scenarios, available in a variety of formats, including Bioclim variables. These are derived from 
the monthly temperature and rainfall values in order to generate biologically meaningful variables, which 
are often used in species distribution modeling and related ecological modeling techniques. The bioclimatic 
variables (Bio01-Biol9) represent annual trends (e. g., mean annual temperature, annual precipitation) 
seasonality (e. g., annual range in temperature and precipitation) and extreme or limiting environmental 
factors (e. g., temperature of the coldest and warmest month, and precipitation of the wet and dry quarters). 
For future climate assessments we used the climate model of 2030 generated for the A1B scenario for emissions 
of greenhouse gases and sulphate aerosols. For tracking retrospective climate change the MERRAclim global set 
of satellite-based bioclimatic variable that have been built for each of the last three decades (1980s, 1990s and 
2000s) was employed (Vega et al., 2017). Grids were presented in a 10 arcmin resolution. 

Co-linearity between predictors may lead to overfitting and impeded predictive performance (Guisan, 
Thuiller, 2005), therefore variable selection analyses was employed by retaining only variables with a Spearman’s 
pair-wise rank correlation coefficient rho < |0.9|. When two variables were highly correlated we chose the one 
least correlated to others, leading to a total of thirteen selected climatic variables (table 1). 


Table 1. Selected climatic variables 














Variable Code Variable 

Bio02 Mean diurnal temperature range (mean(period max-min)) (°C) 

Bio03 Isothermality (Bio02+Bio07) 

Bio05 Max temperature of warmest week (°C) 

Bio06 Min temperature of coldest week (°C) 

Bio07 Temperature annual range (Bio05-Bio06) (°C) 

Bio08 Mean temperature of wettest quarter (°C) 

Bio09 Mean temperature of driest quarter (°C) 

Biol2 Annual precipitation (mm) 

Biol4 Precipitation of driest week (mm) 

Biol5 Precipitation seasonality (C of V) 

Biol7 Precipitation of driest quarter (mm) 

Biol18 Precipitation of warmest quarter (mm) 

Biol9 Precipitation of coldest quarter (mm) 
Results 


The final model was selected based on the fine tuning using the AICc value and the best 
model only used hinge features with a regularization multiplier value of 2.0. The performance 
of the Maxent model is usually evaluated by the threshold-independent receiver operating 
characteristic (ROC) approach, where the calculated area under the ROC curve (AUC) is 
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Fig. 2. Monthly precipitation (mm) change throughout the year. 


considered as a measure of prediction success. The ROC curve is a graphical method that 
represents the relationship between the false-positive fraction (one minus the specificity) 
and the sensitivity for a range of thresholds. It has a range of 0-1, with a value greater than 
0.5 indicating a better-than-random performance event (Fielding, Bell, 1997). A rough 
classification guide is the traditional academic point system (Swets, 1988): poor (0.5-0.6), 
fair (0.6-0.7), good (0.7-0.8), very good (0.8-0.9) and excellent (0.9-1.0). The ROC curve 
of the final model shows a ‘good’ predictive power with an average training AUC score of 
0.733 and standard deviation of 0.023. 

Based on permutation importance, Bio09 (Mean temperature of driest quarter) was 
the most significant variable (36.0 %) followed by Biol4 (Precipitation of driest week) 
contributing 24.1 %; the contribution of Bio12 (Annual precipitation), 9.4 %, was close to 
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Fig. 3. Monthly temperature (°C) change throughout the year. 


476 V. M. Tytar, I. I. Kozinenko, S. V. Mezhzherin 





0.55} 


0.5} 


HS 


0.45} 


0.4} 











-0.24 -012 0 0.12 0.24 036 048 06 0.72 
Mean temperature of driest quarter 
Fig. 4. The marginal response curve for the explanatory variable Bio09 (Mean temperature of driest quarter). 
(HS — habitat suitability). 


the accepted threshold. Among the thirteen selected environmental variables, the first two 
mentioned above accounted for around 60 % of the model prediction. 

The results of the Maxent model’s internal jackknife test of factor importance were 
mostly consistent with the permutation importance. The environmental variable with 
highest gain when used in isolation is Bio12 (Annual precipitation), which therefore appears 
to have the most useful information by itself, and the environmental variable that decreases 
the gain the most when it is omitted is Bio09 (Mean temperature of driest quarter), which 
therefore appears to have the most information that isn’t present in the other variables. 

In terms of the bioclimatic niche of the Steppe mouse two of the emphasized above 
variables seem to play a key role in shaping the ecology and distribution of the species: 











42 48 54 6 66 72 78 84 9 96 
Precipitation of driest week 


Fig. 5. The marginal response curve for the explanatory variable Biol14 (Precipitation of driest week). (HS — 
habitat suitability). 
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Fig. 6. A current climate habitat suitability map for the Steppe mouse (Mus spicilegus) in Ukraine. Darker shades 
of gray denote areas of higher predicted habitat suitability probabilities (= 0.5) and lighter shades correspond 
to lower (> 0.311 and < 0.5). [Administrative regions in Ukraine: 1 — Chernihiv Region; 2 — Kyiv Region; 3 — 
Ternopil Region; 4 — Ivano-Frankivsk Region.] 


namely Bio09 (Mean temperature of driest quarter) and Biol 4 (Precipitation of driest week). 
Both in conjunction can be considered as proxies for describing wintering conditions as far 
as the driest season in the study area falls on January-March when temperatures by large 
exhibit their lowest levels (figs 2 and 3, respectively; the graph in fig. 2 has been smoothed 
using the three-point moving average). 

Habitat preference is interpreted from the response of the predicted habitat suitability 
to a marginal change in each variable, all other variables set to their average value. Values 
giving a high habitat suitability (= 0.5) can be considered as preferable for the species, 
although they may also be found in areas with lower habitat suitability (Wangen et al., 
2016). The response curves for Bio09 (Mean temperature of driest quarter) and Biol4 
(Precipitation of driest week) are presented in fig. 4 and fig. 5, respectively. Because Bio12 
(Annual precipitation) and Biol4 (Precipitation of driest week) are loosely correlated 
(Spearman's rho = 0.76, p = 0.00) their effects on habitat suitability appear in a similar 
fashion, therefore the graph for the former variable is not shown. The y-axis baseline in the 
graphs is set to the 10th percentile training presence logistic threshold value of 0.3708. 

The response to Bio09 (Mean temperature of driest quarter) shows a threshold 
increase from low to high habitat suitability as the mean temperature increases from 
-0.32 °C and reaches a plateau of 0.57 at the point of 0.78 °C; high habitat suitability 
(= 0.5) can be considered to occur in areas where temperatures are above -0.08 °C. On the 
contrary, Bio14 (Precipitation of driest week) demonstrates a downward and close to linear 
trend, suggesting that increasingly wetter conditions lead to lower habitat suitability. Using 
a linear regression model for approximation, we can assume that habitat suitability will 
reach the 10th percentile training presence logistic threshold at around the point of 9.9 mm 
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of precipitation; high habitat suitability (= 0.5) can be considered to occur in areas where 
there is 7.7 or less millimeters of precipitation/driest week in the year. 


Current Distribution 

The probability maps obtained from the 25 replicated runs were averaged to obtain a 
habitat suitability map for the Steppe mouse corresponding to current climate (fig. 6). The 
choice of the cut-off really depends upon the purpose of modelling exercise. For creating 
a range map the threshold defined in the minimum training presence statistics ensures 
that the model captures all training data. This is effectively a species potential distribution: 
everywhere the species could survive in the environment based on the samples (Jobe, Zank, 
2008). For the mean of the 25 minimum training presence thresholds we computed lower 
and upper limits for 95 % confidence intervals, using 9999 bootstrap replicates, and picked 
from there a value of that best suited our purpose in terms of capturing the training data. In 
this respect a threshold value of 0.311 captured 98.6 % of our data. 


Discussion 


In this study we made an attempt to highlight important variables of the current 
bioclimatic niche of the Steppe mouse, realizing that the persistence of populations of 
the species is affected by a combination of physical, biological factors, and anthropogenic 
drivers, which interact in a complex manner. Nevertheless, our coarse scale modeling 
exercise has shown the importance and significance of bioclimatic variables, particularly 
those characterizing wintering conditions, in shaping the niche of the species, although 
model performance (based on the AUC) reaches at the best only ‘good’ predictive power. 
However, such models are considered useful and acceptable for the purpose they are 
intended for (Rykiel, 1996). This occurs when ecological parameters are omitted from the 
modelling framework and lead to the insufficient description of the species’ distribution 
and niche (Hanspach et al., 2010). In our case this was expected because intentionally only 
bioclimatic parameters were used. Of course, it is reasonable to suggest that factors other 
than climate shape the distributions and niche of the species. 

Climatic factors, which normally underlie population change, are expected to exert 
their influence most strongly during certain critical for the survival of the species seasons 
(e. g., the summer dry period or the winter cold period). Over-winter survival is a crucial 
factor in temperate regions, especially for small mammals (Aars, Ims, 2002). Among the 
ecological factors, low ambient temperatures seem to be especially important (Jackson 
et al., 2001). Unlike house mice, mound-building mice live in agricultural fields and are 
found in the steppe zone in Central and Eastern Europe, which have extremely cold but dry 
winters (Simeonovska-Nikolova, 2000). Consequently temperature should be an important 
factor, and the modelling exercise confirms this. In this respect it can be suggested that 
mounds have an insulating role and indeed they reduce temperature variation of the soil 
(Szenczi et al., 2011) 

On the other hand, winters with mild weather seem to be particularly detrimental 
by causing flooding during periods of spring thaw (Aars, Ims, 2002). Similar critical 
weather episodes in winter and early spring have also been suggested to dramatically 
reduce winter survival in small mammals. For instance, deer mouse (Peromyscus 
maniculatus) populations in Colorado responded favourably to rainfall during warm 
periods, but crashed when high rainfall occurred during cold periods (Calisher et al., 
2005; Luis et al., 2010). In our case winter survival of the species most likely depends 
on the combination of temperature and precipitation, when the ambient temperature 
is not too low (to avoid freezing), but not too high (to perhaps trigger snow melting), 
and less precipitation the better. This echoes with A. I. Muntyanu’s view on main 
factors responsible for mouse deaths: low temperature and frequent thaws, leading to 
drenched food and nest chambers (Muntyanu, 1990); from the ENM we can now say 
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Fig. 7. 0.5 °C isotherms for Bio09 (Mean temperature of driest quarter) for different time periods: 1 — 1980s; 
2 — 2000s; 3 — contemporary; 4 — predicted for 2030. 


that this is most likely to happen if there are provisions for a considerable source of 
melt water and/or rain. 

Overwintering structures of the Steppe mouse (mounds) are built from soil and a 
considerable amount of plant material. Recent studies presume that the stored vegetable 
matter is not, or not exclusively, serving as food and indicate that the mounds have not 
only an insulating role, but also waterproofing properties (Hölzl et al., 2011; Szenczi et al., 
2011). On this item we tend to agree with the authors, because once winter precipitation 
is so important it seems reasonable to suggest the development of such waterproofing 
properties. In this case the stored vegetable matter could function as a “thatched roof” and 
help to keep the mound and animals inside in a dry condition. 

Interestingly, areas of M. spicilegus habitat suitability are predicted to occur in the north 
beyond the known home range limits of the species, namely in Kyiv and Chernihiv Regions 
(marked 1 and 2 in fig. 6), from where up to now no presences have been recorded. In terms 
of inventory the model presented here may be useful for guiding search for new populations 
of the Steppe mouse and identifying candidate areas for ground validation of the model. 
Similarly, areas of habitat suitability are predicted to occur in the west, also beyond the known 
home range limits of the species. These are located in Ternopil and Ivano-Frankivsk Regions 
(3 and 4in fig. 6) and appear to occur alongside the Dnister River, which may play a role ofan 
ecological corridor for the spread of the species in a westward direction. 

Although factors other than climate can be responsible for shaping the distributions 
of the considered species, most likely it is climate change that has triggered the recent 
expansion of the home range of the Steppe mouse in Ukraine and neighbouring Russia. If 
so, this expansion process will continue together with the changing climate and new records 
of locations of the species may be used for monitoring such change. To describe the pace 
of climate change that the species should track we mapped for different time periods the 
0.5 °C isotherm for Bio09 (Mean temperature of driest quarter), which was distinguished 
as the most influential variable forming the bioclimatic niche of the Steppe mouse (fig. 7). 
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Since the 1980s there has been a significant shift towards the north, which by and large 
could have favoured the expansion of the species. 


Conclusion 


Climatic conditions have been suggested to be the primary factors that limit the 
range expansion of M. spicilegus in Eastern Europe. Using an ENM approach, we 
have shown in this respect the importance and significance of a number of bioclimatic 
variables, particularly those characterizing wintering conditions, under which higher mean 
temperatures enhance habitat suitability, whereas more precipitation leads to an opposite 
effect. The broadly northwards shift of the home range of the species in Ukraine could 
generally be due to the increasing (since the 1980s) mean temperature of the winter season. 
We expect this expansion process will continue together with the changing climate and 
new records of locations of the species may be used for monitoring such change. 
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